# fracture size and center determination, with fracture parameters acquired in code1.py as input
import matplotlib.pyplot as plt
import numpy as np
import math
import random
import copy


def dist(t1, t2):
    dis = math.sqrt((np.power((t1[0] - t2[0]), 2) + np.power((t1[1] - t2[1]), 2) + np.power((t1[2] - t2[2]), 2)))
    return dis


def loadDataSet(fileName, splitChar):
    dataSet = []
    with open(fileName) as fr:
        for line in fr.readlines():
            curline = line.strip().split(splitChar)
            fltline = list(map(float, curline))
            dataSet.append(fltline)
    return dataSet


def proj(center, occurrence):
    if occurrence[1] == 0:
        a = b = 0
    else:
        ab = math.tan(math.radians(occurrence[1]))
        if ab == 0:
            a = b = 0
        elif occurrence[0] == 0:
            a = 0
            b = ab
        elif occurrence[0] == 180:
            a = 0
            b = - ab
        else:
            if 0 < occurrence[0] < 180:
                a = np.sqrt(ab / (1 + (math.tan(math.radians(90 - occurrence[0]))) ** 2))
            else:
                a = - np.sqrt(ab / (1 + (math.tan(math.radians(90 - occurrence[0]))) ** 2))
            b = a * math.tan(math.radians(90 - occurrence[0]))
    d = - a * center[0] - b * center[1] - center[2]
    return a, b, d


def distances(point, a, b, d):  # distance and project point calculation
    distance = abs(a * point[0] + b * point[1] + point[2] + d) / np.sqrt(a ** 2 + b ** 2 + 1)
    t = - (a * point[0] + b * point[1] + point[2] + d) / (a ** 2 + b ** 2 + 1)
    propit = [point[0] + a * t, point[1] + b * t, point[2] + t]
    return distance, propit


def dbscan(Data, Eps, MinPts):
    """
    DBSCAN for both micro-seismic denoising and fracture size determination
    :param Data: points for denoising or clustering
    :param Eps:
    :param MinPts:
    :return:
    """
    num = len(Data)
    unvisited = [i for i in range(num)]
    visited = []  # 已经访问的点的列表
    C = [-1 for _ in range(num)]
    k = -1
    cluster, tot_cluster, reclus = [], [], []
    cluspoints = 0
    while len(unvisited) > 0:
        subcluster = []
        p = random.choice(unvisited)
        unvisited.remove(p)
        visited.append(p)
        N = []
        for i in range(num):
            if dist(Data[i], Data[p]) <= Eps:
                N.append(i)
        if len(N) >= MinPts:
            k += 1
            C[p] = k
            for pi in N:
                if pi in unvisited:
                    unvisited.remove(pi)
                    visited.append(pi)
                    M = []
                    for j in range(num):
                        if dist(Data[j], Data[pi]) <= Eps:
                            M.append(j)
                    if len(M) >= MinPts:
                        for t in M:
                            if t not in N:
                                N.append(t)
                if C[pi] == -1:
                    C[pi] = k
            gt = 0
            for jc in N:
                if Data[jc] not in cluster:
                    cluster.append(Data[jc])
                    gt += 1
                    if Data[jc] not in subcluster:
                        subcluster.append(Data[jc])
            cluspoints += gt
            tot_cluster.append(subcluster)
        else:
            C[p] = -1
            if Data[p] not in reclus:
                reclus.append(Data[p])
    clustered = []  # only one cluster can exist in a fracture, so the largest cluster is selected here
    for jt in range(len(tot_cluster)):
        if jt == 0:
            clustered = tot_cluster[jt]
            ma = len(tot_cluster[jt])
            continue
        elif len(tot_cluster[jt]) > ma:
            ma = len(tot_cluster[jt])
            passed = clustered
            clustered = tot_cluster[jt]
        else:
            passed = tot_cluster[jt]
        if passed:
            for ct in passed:
                if ct not in reclus:
                    reclus.append(ct)
    if not len(clustered) + len(reclus) == len(Data):
        for jtc in clustered:
            if jtc in reclus:
                reclus.remove(jtc)
    return cluster, clustered, reclus


def genfrasize(sigma, avery, avery2, points):
    """
    generating fracture size from micro-seismic events and fracture parameters
    :param sigma: distance threshold
    :param avery: search radius for DBSCAN micro-seismic denoise
    :param avery2: search radius for DBSCAN fracture size determination (attempts should be made to determine)
    :param points: minimum number of points (defines the same in the both DBSCAN methods)
    :return: fracture fit result (tfit), unfitted points (unf)
    """
    clcc = []
    with open('list.txt') as fr:
        for line in fr.readlines():
            curline = line.strip().split(',')
            clcc.append(curline[0])
    rawdataa = loadDataSet(f'{clcc[0]}.txt', ',')
    raw = []
    for ab in rawdataa:
        if ab not in raw:
            raw.append(ab)
    data, sudata, rec = dbscan(raw, avery, points)
    paras = loadDataSet('best_fit.txt', ',')
    para = []
    for jt in paras:  # convert fracture parameters to analytical form
        a, b, d = proj(jt[:3], jt[3:])
        spara = [a, b, d]
        para.append(spara)
    redata = copy.deepcopy(data)
    fitp, fitd = [[] for i in range(len(redata))], [[] for i in range(len(redata))]
    dsts = [10 ** 15 for ip in range(len(redata))]
    tpts = []
    for ij in para:
        for tj in redata:
            dst, gec = distances(tj, ij[0], ij[1], ij[2])
            fitd[redata.index(tj)].append(dst)
            if para.index(ij) == 0:
                tpts.append(tj)
            if para.index(ij) == 0 or dst < dsts[redata.index(tj)]:
                dsts[redata.index(tj)] = dst
    unf, sfitd, sfitp = [], [], []
    for cac in redata:
        if dsts[redata.index(cac)] > sigma:  # select unfitted points
            unf.append(cac)
            fitd.remove(fitd[tpts.index(cac)])
            tpts.remove(cac)
    print(f'{len(tpts)} out of {len(redata)} points fitted, fitted rate={len(tpts) / len(redata) * 100}%')
    tfit = []
    for dcc in fitd:
        dccs = sorted(dcc)  # sort the distances value from small to large
        ssfitp = []
        for scc in dccs:
            ssfitp.append(dcc.index(scc) + 1)  # fracture number should also be changed accordingly
        sfitd.append(dccs)
        sfitp.append(ssfitp)
    btpts = copy.deepcopy(tpts)
    tofit, rejected = [], []
    numrj = [0 for j in range(len(tpts))]  # rejects for each point
    for a in range(len(para)):
        afit = []
        for b in range(len(tpts)):
            if sfitp[b][0] == a + 1:
                afit.append(btpts[b])
        tofit.append(afit)  # temporary classification result
    for afa in tofit:
        af, cafit, rej = dbscan(afa, avery2, points)
        rejected.append(rej)  # "no home" points that are to be classified, classified according to fracture
        tfit.append(cafit)
    brjt, cfra = [], []
    for dcd in rejected:
        for ede in dcd:
            brjt.append(ede)
    for bcb in tpts:  # merge the "no home" points
        if bcb in brjt:
            cfra.append(sfitp[tpts.index(bcb)])
            numrj[tpts.index(bcb)] += 1
    print(f'{len(brjt)} out of {len(redata) - len(unf)} points remaining, max rejections={max(numrj)}')
    # repeat above step until all the points are classified or the number of attempts reaches to 4 times of the
    # number of fractures
    while len(brjt) > 0 and max(numrj) < 4 * len(para):
        rejected, tofit = [], []
        for aca in range(len(tfit)):
            for bdb in brjt:
                if cfra[brjt.index(bdb)][numrj[tpts.index(bdb)] % len(cfra[brjt.index(bdb)])] == aca + 1:
                    tfit[aca].append(bdb)
        for afas in tfit:
            af, cafit, rej = dbscan(afas, avery2, points)
            rejected.append(rej)
            tofit.append(cafit)
        brjt, cfra = [], []
        for dcd in rejected:
            for ede in dcd:
                brjt.append(ede)
        for bcb in brjt:
            cfra.append(sfitp[tpts.index(bcb)])
            numrj[tpts.index(bcb)] += 1
        tfit = copy.deepcopy(tofit)
        print(f'{len(brjt)} out of {len(redata) - len(unf)} points remaining, max rejections={max(numrj)}')
    return tfit, para, unf, redata, len(brjt)


def defrasize(tfit, para, unf, data, cva):
    """
    fracture size, center and vertices' calculation from classification results
    :param tfit: classification results
    :param para: analytical parameters of fractures
    :param unf: unfitted points
    :param data: effective micro-seismic events
    :param cva: number of "no home" points
    :return:
    """
    fig = plt.figure(figsize=(8, 6), dpi=80)
    ax = fig.add_subplot(111, projection='3d')
    file3 = open('define_4_user_ellipses.dat', 'a')
    # write input according to the format of the corresponding LaGriT input file
    file3.write(f'nEllipses: {len(tfit)}' + '\n')
    file3.write('nNodes: 4' + '\n')
    file3.write('Coordinates:' + '\n')
    fits = 0
    for t in range(len(para)):
        real_fitted = tfit[t]
        x, y, z = [], [], []
        x_tot = y_tot = f_1 = f_2 = 0
        for tjs in real_fitted:
            cdc, prjp = distances(tjs, para[t][0], para[t][1], para[t][2])
            x.append(prjp[0])
            y.append(prjp[1])
            z.append(prjp[2])
            x_tot += prjp[0]
            y_tot += prjp[1]
            f_1 += prjp[0] * prjp[1]
            f_2 += prjp[0] ** 2
        x_tot /= len(x)
        y_tot /= len(y)
        # least square method for fitting central line
        slope = (f_1 - len(x) * x_tot * y_tot) / (f_2 - len(x) * x_tot ** 2)
        # degree of the line
        alpha = math.degrees(math.atan(slope))
        if slope < 0:
            alpha += 180
        # rotation angle determination
        if 0 <= slope <= 1:
            theta = alpha
        elif -1 <= slope < 0:
            theta = alpha - 180
        else:
            theta = alpha - 90
        xj, yj = [], []
        # coordinate transformation
        for j in range(len(x)):
            x_new = (x[j] - x_tot) * math.cos(math.radians(theta)) + (y[j] - y_tot) * math.sin(math.radians(theta))
            y_new = (y[j] - y_tot) * math.cos(math.radians(theta)) - (x[j] - x_tot) * math.sin(math.radians(theta))
            xj.append(x_new)
            yj.append(y_new)
        x, y, z = [], [], []
        for tjs in real_fitted:
            x.append(tjs[0])
            y.append(tjs[1])
            z.append(tjs[2])
        ax.scatter(x, y, z, marker='o', s=60, alpha=1., label='microseismic events')
        xt, yt = np.linspace(min(xj), max(xj), 15), np.linspace(min(yj), max(yj), 15)
        Xt, Yt = np.meshgrid(xt, yt)
        # inverse coordinate re-transformation
        Ztx = Xt * math.cos(math.radians(theta)) - Yt * math.sin(math.radians(theta)) + x_tot
        Zty = Xt * math.sin(math.radians(theta)) + Yt * math.cos(math.radians(theta)) + y_tot
        Ztz = - para[t][0] * Ztx - para[t][1] * Zty - para[t][2]
        # the 4 vertices of a single fracture
        v_x = [Ztx[0][0], Ztx[0][len(Ztx[0]) - 1], Ztx[len(Ztx) - 1][0], Ztx[len(Ztx) - 1][len(Ztx[len(Ztx) - 1]) - 1]]
        v_y = [Zty[0][0], Zty[0][len(Zty[0]) - 1], Zty[len(Zty) - 1][0], Zty[len(Zty) - 1][len(Zty[len(Zty) - 1]) - 1]]
        v_z = []
        xcot = ycot = zcot = 0
        for i in range(len(v_x)):
            zzz = - para[t][0] * v_x[i] - para[t][1] * v_y[i] - para[t][2]
            zcot += zzz
            v_z.append(zzz)
            xcot += v_x[i]
            ycot += v_y[i]
        center = [xcot / 4, ycot / 4, zcot / 4]  # center of fracture
        vector_2 = [v_x[1] - v_x[3], v_y[1] - v_y[3], v_z[1] - v_z[3]]
        lv_2 = np.sqrt(vector_2[0] ** 2 + vector_2[1] ** 2 + vector_2[2] ** 2)
        vector_3 = [v_x[1] - v_x[0], v_y[1] - v_y[0], v_z[1] - v_z[0]]
        lv_3 = np.sqrt(vector_3[0] ** 2 + vector_3[1] ** 2 + vector_3[2] ** 2)
        scalar = vector_3[0] * vector_2[0] + vector_3[1] * vector_2[1] + vector_3[2] * vector_2[2]
        eqradius = np.sqrt(np.sqrt(lv_3 ** 2 * lv_2 ** 2 - scalar ** 2) / math.pi)  # equilibrium radius
        ax.plot_surface(Ztx, Zty, Ztz, rstride=1, cstride=1, alpha=.7, label=f'fracture {t + 1}')
        print(f'fracture {t + 1} is created, equal radius={eqradius}, center={center}, '
              f'with {len(real_fitted)} points fitted')
        # write the 4 vertices of a fracture to file
        file3.write('{' + f'{v_x[0]}' + ',' + f'{v_y[0]}' + ',' + f'{v_z[0]}' + '}' + '\t')
        file3.write('{' + f'{v_x[1]}' + ',' + f'{v_y[1]}' + ',' + f'{v_z[1]}' + '}' + '\t')
        file3.write('{' + f'{v_x[3]}' + ',' + f'{v_y[3]}' + ',' + f'{v_z[3]}' + '}' + '\t')
        file3.write('{' + f'{v_x[2]}' + ',' + f'{v_y[2]}' + ',' + f'{v_z[2]}' + '}' + '\n')
        fits += len(real_fitted)
    file3.close()
    x, y, z = [], [], []
    for tjt in unf:
        x.append(tjt[0])
        y.append(tjt[1])
        z.append(tjt[2])
    print(f'fitted rate={100 - len(x) / len(data) * 100}%')
    print(f'real draw={fits / (len(data) - cva) * 100}%')
    ax.scatter(x, y, z, c='k', marker='o', s=10, alpha=1., label='unfitted events')
    ax.set_xlabel('X')
    ax.set_ylabel('Y(N)')
    ax.set_zlabel('Z')
    # coordinate and view angle options (may differ because of differences in micro-seismic events and fractures)
    ax.set_zlim(-4400, -3200)
    ax.set_xlim(-600, 600)
    ax.set_ylim(-600, 600)
    ax.view_init(azim=-66, elev=20)
    plt.show()
    return


def approximate465(anyj, t):
    if anyj * 10 ** t - int(anyj * 10 ** t) > 0.5:
        anyja = (int(anyj * 10 ** t) + 1) / 10 ** t
    elif anyj - int(anyj) < 0.5:
        anyja = (int(anyj * 10 ** t)) / 10 ** t
    else:
        if (int(anyj * 10 ** t) + 1) % 2 == 0:
            anyja = (int(anyj * 10 ** t) + 1) / 10 ** t
        else:
            anyja = (int(anyj * 10 ** t)) / 10 ** t
    return anyja


opt = loadDataSet('option.txt', ',')
sigma, avery, points, iters = opt[1][0], opt[1][1], opt[1][2], int(opt[1][3])
avery2 = eval(input('please input a threshold for size determination:'))  # input DBSCAN parameter
tfit, para, unf, ini, ctb = genfrasize(sigma, avery, avery2, points)
defrasize(tfit, para, unf, ini, ctb)
